% 上图: LE1, LE2, LE3 随 mu 的曲线
% 下图: x 的分岔图

clc; clear;

k = 0.35;                            % 固定 k
mu_values = linspace(2.6, 3.6, 350); % mu 扫描范围，采样 1000 个点
x0 = 0.1; y0 = 0; z0 = -1;           % 初始条件
nTransient = 500;                    % 丢弃瞬态步数
nLE = 100;                           % 用于计算 LE 的步数
nObserve = 500;                      % 分岔图观测步数
nStore = 100;                        % 分岔图中每个 mu 存储多少点

% 预分配存储数组
% 用于存储 LE
LE1 = zeros(size(mu_values));
LE2 = zeros(size(mu_values));
LE3 = zeros(size(mu_values));
% 用于存储分叉
bif_mu = zeros(length(mu_values)*nStore, 1);
bif_x  = zeros(length(mu_values)*nStore, 1);
idx = 1;

% 扫描mu
for i = 1:length(mu_values)
    mu = mu_values(i);

    
    state = [x0, y0, z0];

    % (a) 丢弃瞬态
    for n = 1:nTransient
        [dx, dy, dz] = mclm(state, mu, k);
        state = [dx, dy, dz];
    end

    % (b) 计算 Lyapunov 指数
    LE_vals = LEs(state, mu, k, nLE);  % 返回 [LE1, LE2, LE3]
    LE1(i) = LE_vals(1);
    LE2(i) = LE_vals(2);
    LE3(i) = LE_vals(3);

    % (c) 分岔图：再做 nStore 次迭代，记录 x
    for n = 1:nObserve
        [dx, dy, dz] = mclm(state, mu, k);
        state = [dx, dy, dz];
        
        % 存储最后nStore个点
        if n > (nObserve - nStore)
            bif_mu(idx) = mu;
            bif_x(idx)  = dx;  % 使用更新后的dx（即x_{n+1}）
            idx = idx + 1;
        end
    end
end


% 上方：LE1, LE2, LE3 随 mu 的变化
subplot(2,1,1);
plot(mu_values, LE1, 'r', 'LineWidth',1.2); hold on;
plot(mu_values, LE2, 'g', 'LineWidth',1.2);
plot(mu_values, LE3, 'b', 'LineWidth',1.2);
xlabel('\mu');
ylabel('LE');
legend('LE_1','LE_2','LE_3','Location','best');
title(sprintf('LE Spectra'));
grid on;

% 下方：分岔图 (bifurcation diagram)
subplot(2,1,2);
scatter(bif_mu, bif_x, 1, 'r', 'filled');
xlabel('\mu'); ylabel('x');
title('Bifurcation');
grid on;
